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Multiple mass solvers 

B. Jegerlehner^ 

'^Indiana University, Depatrmcnt of Physics, Swain Hall West 117, 
Bloomington, IN 47408, USA 

We present a general method to construct multiple mass solvers from standard algorithms. As an example, the 
BiCGstab-M algorithm is derived. 



1. INTRODUCTION 

It has been discussed recently that, using 
Krylov space solvers, the solutions of shifted lin- 
ear equations, where {A + a)x ~ b — has to 
be calculated for a whole set of values of a, can 
be found at the cost of only one inversion. This 
kind of problem arises in quark propagator cal- 
culations for QCD as well as other parts of com- 
putational physics (see It has been realized 
that several algorithms allow to perform this task 
using only as many matrix-vector operations as 
the solution of the most difficult single system re- 
quires. This has been achieved for the QMR 0, 
the MR 1^ and the Lanczos-implementation of 
the BiCG method ^ . We present here a unifying 
discussion of the principles to construct such al- 
gorithms and succeed in constructing shifted ver- 
sions of the CG, CR, BiCG and BiGGstab algo- 
rithms [Q, using only two additional vectors for 
each mass value. 

2. METHOD 

The iterates of Krylov space methods, espe- 
cially the residuals r„, are generally polynomials 
P„ (A) of the matrix applied to some initial vector 
Tq. The polynomial generating the residual gen- 
erally has to be an approximation to zero in some 
region enclosing the spectrum of A while statis- 
fying P„(0) = 1. The key to the method is the 
observation that shifted polynomials, defined by 



P,:(A + a)=<P„(A), 



(1) 



multiple (T values using no additional matrix- 
vector multiplications. We expect cj^ < 1 if the 
condition number of ^ + ct is smaller than the one 
of A, which is confirmed in numerical tests. 

Generally the polynomial generated in a solver 
is defined by some recursion relation. We will 
therefore need to know the recursion relation for 
the shifted polynomial, too, which can be found 
easily by parameter matching. Here, we discuss 
only polynomials which statisfy P(0) = 1, but 
more general normalisation conditions are han- 
dled analogously. For the two-term recursion re- 
lation 



P„+i(.t) = {anX + l)Pn{x) 

wc find for the polynomial shifted by a 

-fri+l(^) 




are useful objects, since vectors generated by 
these shifted polynomials can be calculated for 



(2) 

(3) 
(4) 



This formula has also been found in (2| with dif- 
ferent methods. From ^ we can read off the 
parameters of the shifted polynomial, while (jj) 
determines cj^. Note that this formula holds for 
any choice of «„. We can easily generalize this 
method to three-term recurrences 

P,i+l{x) = {anX + (3n)Pn{x) + 7„P„-l(x) (5) 

and find explicit expressions for the parameters 
of the shifted polynomial. It turns out that the 
Lanczos polynomials for the matrices A and A+a 
fulfill (|l|) automatically (this was the original ob- 
servation in We derived the shifted polyno- 
mial however for arbitrary choices of parameters 



2 



a and /3. The most interesting case are coupled 
two-term recurrences, since they have superior 
stability properties over three-term recurrences. 
We consider the general recurrence 



Qn{x) 
Pn+l{x) 



(6) 
(7) 



In CG-type algorithms, the parameters are cho- 
sen so that Pn{x) is the Lanzcos-polynomial (nor- 
mahzed to P(0) = 1). We thus demand P.^{x + 
a) = d^Pnix). By transforming the above rela- 
tion to a simple three-term recurrence and apply- 
ing the formulae found in this case for the shifted 
parameters we find 



+ (1 - fT/3„) 



(8) 

(9) 
(10) 



It can easily be checked that Q'^^{x + a) ^ cQn{x), 
so if we want to use this recursion relation in an 
algorithm we have to replace 



f3:^{x + a)QUx) 



iF„+i(x)-<P„(a;) (11) 



in the shifted systems. Since in CG-type algo- 
rithms the update of the solution vector involves 
Qn{x)vQ, this vector has to be iterated and stored 
for all shifted systems. 

3. APPLICATIONS 

Using the above formulae, we can easily de- 
rive a variety of linear system solvers as shown 
in table |l|. We present here only the algorithm 
of greatest interest for quark propagator calcula- 
tions, BiCGstab-M. 

3.1. BiCGstab-M 

The BiCGstab-M algorithm is a mixture be- 
tween the BiCG and the MR algorithm. It is 
therefore not surprising that we can simply use 
the formulae for the two-term and the coupled 
two-term recurrences and construct a shifted al- 
gorithm. In the BiCGstab algorithm we gen- 



Method 



MR-M 
CR-M 
QMR-M (3-term) 
QMR-M (2-term) 
TFQMR-M 
BIORESU 
BiCG-M 
BiCGstab-M 



Referc 
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|6 



111 



|3 



|41 



4] 



4] 



Memory 



N 
2N 

m 

3N 
5N 
2N 
2N 
2N+1 



Table 1 

Memory requirements and references for shifted 
system algorithms for unsymmetric or nonher- 
mitean matrices. We list the number of addi- 
tional vectors neccessary for N additional values 
of (7 ( which is independent of the use of the 75- 
symmetry) . 



erate the following sequences 
r„ - P„(A)i?„(A)ro (12) 

Wn = Pn{A)Rn-l{A)ro (13) 

Sn = Qn{A)Rn{A)r„ (14) 

where Zn{x) and Qn{x) are exactly the BiCG- 
polynomials and _R„ (x) is derived from a minimal 
residual condition. For the shifted algorithm we 
demand 

P^ix + a) = c^P„(x), Kix + <J) = d^i?„(a;).(15) 

Using the above formulae we can explicitly de- 
termine the constants c and d and the shifted 
parameters ot the polynomials. The remaining 
difficulty is to derive the iteration for the solu- 
tion Xn and the vector s„. The update of these 
two vectors has the form 

Xn+l = Xn~ (3nS„ + XnWn+1 (16) 
Sn+1 = rn+1 + an+liSn - XnASn) (17) 

This means we have to eliminate Asn from the 
update of s„. The updates for the shifted vectors 
x'^ and sj^ then look as follows: 

xZ~l3'^X, + x^cy^,_,w^+i 



(18) 
(19) 



t-'n 
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We therefore need to introduce 2 vectors for each 
shifted system and one additional vector to store 
r„. Note that the case /?„ = leads to a break- 
down of the BiCGstab algorithm and does not in- 
troduce any new problems for the shifted method. 

The convergence of the shifted algorithms can 
be verified by checking that Cnr„ < 1 . It is how- 
ever generally advisable for all shifted algorithms 
to test all systems for convergence after the algo- 
rithm finishes since a loss of the conditions (15) 
due to roundoff errors might lead to erratic con- 
vergence. 

4. LIMITATIONS 

The most serious limitation of the method is 
given by the fact that the starting residual for 
all systems must be the same, which excludes 
cr-dependent left preconditioning. Furthermore, 
preconditioning must retain the shifted structure 
of the matrix. This means that especially even- 
odd preconditioning is not applicable. To stabi- 
lize the algorithm, however, one can apply poly- 
nomial preconditioning: 



P[A)Ay = 6, P[A)y. 
We must have 

P''{A + <t){A + g)=P{A)A + ti. 



(20) 



(21) 



Note that P'^ might not be a good preconditioner 
for A + a, which is compensated for by the faster 
convergence of the shifted system. A linear pre- 
conditioner, which was also proposed in is 
given by 



(22) 



For the Wilson and clover matrix, this polynomial 
has the property that P'^{x) is a good precondi- 
tioner for A + a. This preconditioner has been 
found to work well in those cases. Higher order 
polynomials can be derived from condition (pT|). 

5. CONCLUSIONS 

We presented a simple point of view to under- 
stand the structure of Krylov space algorithms for 
shifted systems, allowing us to construct shifted 
versions of most short recurrence Krylov space 



algorithms. The shifted CG-M and CR-M al- 
gorithm can be applied to staggered fermion 
calculations. Since efficient preconditioners for 
the staggered fermion matrix are not known, a 
very large improvement by these algorithms can 
be expected. We presented the BiCGstab-M 
method, which, among the shifted algorithms, is 
the method of choice for quark propagator calcu- 
lations using Wilson (and presumably also clover) 
fermions if enough memory is available. The nu- 
merical stability of the algorithms has been found 
to be good [Q . Roundoff errors might however in 
some cases affect the convergence of the shifted 
systems so that the final residuals have to be 
checked. Other discussions can be found in 
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